Isostaticity at Frictional Jamming 
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Amorphous packings of frictionless, spherical particles are isostatic at jamming onset, with the 
number of constraints (contacts) equal to the number of degrees of freedom. Their structural and 
mechanical properties are controlled by the interparticle contact network. In contrast, amorphous 
packings of frictional particles are typically hyperstatic at jamming onset. We perform extensive 
numerical simulations in two dimensions of the geometrical asperity (GA) model for static friction, 
to further investigate the role of isostaticity. In the GA model, interparticle forces are obtained 
by summing up purely repulsive central forces between periodically spaced circular asperities on 
contacting grains. We compare the filling fraction, contact number, mobility distribution, and 
vibrational density of states using the GA model to those generated using the Cundall- Strack (CS) 
approach, where static friction is modeled by a tangential linear spring when two disks form a contact 
that is allowed to slide when their relative tangential displacement exceeds the Coulomb threshold. 
We find that static packings of frictional disks obtained from the GA model are mechanically stable 
and isostatic when we consider interactions between asperities on contacting particles. The crossover 
from frictionless to frictional behavior as a function of the static friction coefficient coincides with 
a change in the type of interparticle contacts and the disappearance of a peak in the density of 
vibrational modes for the GA model. These results emphasize that mesoscale features of the model 
for static friction play an important role in determining the structural and mechanical properties of 
granular packings. 



Over the past fifteen years, intense effort has been de- 
voted to understanding the jamming transition in ather- 
mal systems composed of frictionless spheres with purely 
repulsive contact interactions [TJ-Ql - However, physical 
models of granular media must include static friction [E| . 
Several recent experimental 0, 0] and simulation stud- 
ies [|-[Iq| have provided significant insight into the jam- 
ming transition for frictional, spherical grains. Amor- 
phous, static packings of frictional spheres can be ob- 
tained at jamming onset over a wide range of contact 
numbers, d + 1 < z < 2d 0, 0j], 03, wnere d is the spa- 
tial dimension. In addition, amorphous sphere packings 
show a crossover from frictionless behavior with pack- 
ing fraction near random close packing <\> ~ </>rcp and 
z ~ 2d to frictional behavior with packing fraction near 
random loose packing qb ~ </>rlp and z ~ d + 1 as the 
static friction coefficient increases above a characteristic 
value /i* - 0.1 (0.01) in 2D (3D) [l2|. These prior studies 
have also shown that static packings of frictional spheres 
possess a large number N s of 'sliding' contacts with tan- 
gential forces near the Coulomb yield threshold for small 
/i, and N s decreases with increasing ja [l2|, [T^j. Some 
evidence has been presented that when contact-counting 
arguments account for sliding contacts, frictional sphere 
packings can be described as 'isostatic' with a plateau 
in the density of vibrational modes at low frequencies as 
found for static packings of frictionless spheres fioj ] . 

In this Letter, we address several fundamental open 
questions raised by prior studies of static packings of 
frictional particles: 1) How sensitive are the structural 
and mechanical properties of static packings of frictional 



particles to the friction model employed? 2) What de- 
termines the static friction coefficient /x* that marks the 
crossover from frictionless to frictional behavior for static 
packings? 3) How do the vibrational modes for pack- 
ings of frictional spheres differ from those for static pack- 
ings of frictionless particles with complex and anisotropic 
shapes? 

We employ a mesoscopic, geometrical asperity (GA) 
model for static friction in which the interparticle forces 
are obtained by summing up purely repulsive central 
forces (proportional to overlap) between periodically 
spaced asperities on contacting, otherwise circular disks. 
In contrast, most prior studies have followed the Cundall- 
Strack approach [14[ , where static friction is modeled by 
a tangential spring (with spring constant k t and restor- 
ing force k t u t , where u t is the relative tangential dis- 
placement) when two particles come into contact and the 
Coulomb sliding condition is enforced during the lifetime 
of the contact. The GA model offers several advantages 
for studying static packings. By using the GA model 
we are able to distinguish between interparticle contacts 
based on the number of asperities that interact at each 
contact and calculate the density of vibrational modes 
by taking derivatives of the total potential energy with 
respect to all nontrivial degrees of freedom without mak- 
ing ad hoc assumptions about sliding contacts [10]. Sev- 
eral prior studies have implemented asperity models to 
mimic frictional interactions [HI, [l6|, but these studies 
have mainly focused on dense flows, not static packings. 

Our key finding is that static packings from the GA 
model are mechanically stable and isostatic when we con- 




FIG. 1: Top: Mechanically stable packings of N = 6 bidis- 
perse disks at jamming onset obtained using the CS (left) and 
GA (right) models for friction with /x, /x e ff — 0.3. These pack- 
ings from the CS and GA models are nearly identical with 
cj)j ~ 0.78 and 0.76, respectively. They possess the same 9 in- 
terparticle contacts, and the GA model has the isostatic num- 
ber of contacting asperities N^ a = SN — 1 = 17 for disks in 
2D with three degrees of freedom per disk in periodic bound- 
ary conditions. In the right panel, the central particle has 
five interactions between asperities on three contacting grains. 
The solid and striped gray contacts between the central par- 
ticle and its neighbors represent single and double asperity 
contacts, respectively. Bottom: (left) Schematic of the ratio 
of the tangential and normal forces ft/ fn at constant inter- 
particle overlap versus the relative tangential displacement 
ut for the CS (dashed) and GA (solid) models. For the CS 
model, ft/ f n increases linearly with ut with slope kt, while 

ft I fn — u t / (a if — r^f) 2 — u\ for the GA model, where 

the separation between asperities r^ a is constant at fixed 
overlap. For the GA model, single asperity contacts start 
at ft/ fn = 0, while double asperity contacts give rise to the 
maximal/minimal ft/ fn- The sliding limit for the CS model 
is ±u a = ±/nf n /kt (above which ft/ fn = ±M remains fixed), 
while u a = ±cTij a /(2y/l + 1/Meff) an d ft/ fn is periodic in ut 
with period 4u a in the zero overlap limit for the GA model, 
(right) Schematic of the interaction between two rough disks 
in the GA model with radius R, N a circular asperities with 
radius R a , and angle between bumps 2tt /N a . 



sider interactions between asperities on contacting grains, 
independent of the effective static friction coefficient. We 
also show that the crossover in the structural and me- 
chanical properties of static packings from frictionless to 
frictional behavior as a function of the effective friction 
coefficient coincides with a change in the number of inter- 
actions between asperities per contact and the disappear- 
ance of a strong peak in the density of vibrational modes 
at low frequency with primarily rotational content. Fur- 
ther, we find that the density of vibrational modes for 
the GA model differs from previous calculations for static 
packings of frictional disks using the Cundall-Strack (CS) 
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FIG. 2: Top: Average packing fraction for MS packings 
at jamming onset from the CS and GA models versus the 
static friction coefficient (/i or /i e ff)- The inset in the lower 
left shows ((j>j) versus \i or /i e ff for several system sizes N and 
numbers of asperities N a . The first and second numbers listed 
in the legends correspond to N a and N, and axes without 
tick labels are the same as those in the main panel. Bottom: 
Average interparticle contact number (z pp ) versus the static 
friction coefficient. The insets in the lower left and upper right 
show the N and N a dependence of (z pp ) and the fraction of 
floater particles Nf/N, respectively. 



model [10[. 

We construct mechanically stable packings of N rough 
bidisperse disks (50 — 50 by number with diameter ra- 
tio r = 1.4) in d = 2 using the GA model and com- 
pare their structural and mechanical properties to MS 
packings of frictional disks obtained using the CS ap- 
proach. As shown in the lower right panel of Fig. [1] 
rough circular disks in the GA model are characterized 
by a given number of circular asperities N a with cen- 
ters on the rim of the disks and ratio of the asperity 
to particle radius R a /R- We consider two types of in- 
teractions between disks: 1) an asperity on disk i and 
an asperity on disk j and 2) the core of disk i and an 
asperity on disk j. All interactions have the form of 
purely repulsive linear spring potentials. The interac- 
tion between asperities a and a' on disks i and j is mod- 



eled as Vfi a 
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where r^ a is center-to-center separation between asper- 
ities, aff = Rf + Rf, and = aff + R? + Rf. We 
locate asperity a on the rim of disk i at angle Of = + 
and coordinates = (r^ + cos#f,i*i + sin#f), where 
is the position of disk i. The interaction between 
asperity a on disk i and the core of disk j is mod- 
eled as V* = e/(2af j )(af j - r^-) 2 e(l - r^/a^), where 
cr?- = R? + Ri + Rj . The total potential energy for the 

GA model is V = Ei>j E a >a' + Ei>j E„ V*. 
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FIG. 3: Accumulated mobility distributions for N — 48 bidis- 
perse disks for the CS (left) and GA (right) models for static 
friction coefficients /x, /x e ff = 10" 3 , 10" 1 , 1, and 10, where the 
mobility £ = \ft\/(fJ>fn)- For the GA model, mobilities £ > 1 
can occur due to finite interparticle overlaps. The bin at £ = 1 
includes all £ > 1 to allow a comparison with the CS model. 

For the GA model, we can define an effective static fric- 
tion coefficient as /i e ff = l/- V /((2i? a /i?)/sin(7r/Ar a )) 2 -l, 
which is the ratio of the maximum tangential to nor- 
mal interparticle force when an asperity on disk i fits 
in between two asperities on disk j as shown in the up- 
per right panel of Fig. [H We verified that this is the 
maximum ratio of tangential to normal forces in the 
limit of zero interparticle overlap. We fix the ratio of 
the number of asperities on the large and small parti- 
cles close to r so that the inter-species /i e ff is approx- 
imately the same as that for the large and small par- 
ticles. For the CS model 0, EH, geometrically smooth 
circular disks interact via the purely repulsive linear (or 
Hertzian) spring potential Vij, and static friction is in- 
cluded using a tangential spring with spring constant 
kt/k n = 1/3 (k n = e/aij) [5], and f t remains at its max- 
imum/minimum ±/i/ n when the relative tangential dis- 
placement exceeds the Coulomb threshold. We studied 
system sizes from N = 6 to 96 particles, numbers of as- 
perities N a = 8, 16, and 32, and static friction coefficients 
from /i, /i e ff = 10 -3 to 10. 

For both the GA and CS models, we generate more 
than 10 5 MS packings at jamming onset, for each N 
and /jl or /i e fj, using the compressive-quench-from-zero- 
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FIG. 4: Top: Probability P c of different types of interpar- 
ticle contacts versus the static friction coefficient \i or /x e ff- 
The contact types are single (rightward triangles) and dou- 
ble asperity (leftward triangles) for the GA model and low 
(squares) and high mobility (pentagons) with £ < Q c = 0.5 
and £ > f c , respectively, for the CS model. Bottom: The av- 
erage isostaticity parameter (a) = (((N™ + l)/3 + N f )/N), 
where N£ a is the total number of contacts between asperities, 
versus /x e ff for the GA model for several system sizes N and 
number of asperities N a . (N = 6 and N a = 16, circles; 6 
and 32, squares; 12 and 16, rightward triangles; 12 and 32, 
leftward triangles; 24 and 16, upward triangles; 24 and 32, 
downward triangles; 48 and 16, stars; 48 and 32, hexagons.) 
a = 1 indicates an isostatic number of asperity contacts. 

density simulation protocol [13] • We begin by randomly 
placing disks at zero packing fraction in a square cell of 
unit size with periodic boundary conditions. We increase 
the radius of all particles in small steps corresponding to 
changes of A(f> = 10 -4 . After each packing fraction incre- 
ment, the system is relaxed to the nearest local potential 
energy minimum using dissipative forces proportional to 
the translational and angular velocities of the disks with 
large damping coefficients. If the system after minimiza- 
tion possesses zero total potential energy per particle (i.e. 
less than a small threshold Vt \/e = 10 -14 ), we decom- 
press the system. Otherwise, we compress the system. 
The packing fraction increment is halved each time the 
process switches from compression to decompression or 
vice versa. We perform successive compressions and de- 
compressions until Vtoi < V/N < l.OlVtoh an d the aver- 
age particle overlap is less than 10 _T . For the GA model, 
we find that all packings are mechanically stable with 
Nd = 37V' — 2 eigenvalues > for the dynamical ma- 
trix M ki = d ^ ni , where R = {r x , r 2 , . . . , r N ^ (Ri + 

(R2 + R%)0 2 , - - - , (Rn> + R a N>WN>}, N' = N-N f , 
and Nf is the number of floater particles with fewer than 
three asperity contacts for the GA model. (For the CS 
model, floater particles possess fewer than two interpar- 
ticle contacts.) 

In Fig. [2j we show results for the average packing 
fraction ((f) j) and interparticle contact number (z pp ) = 



3 



(2N PP /(N — Nf)) at jamming onset for the CS and GA 
models. As found in previous studies of static packings 
of frictional disks [12], varies from w 0.84 to 0.75 

and (z pp ) ranges from « 4 to 3 for the CS model as \i 
increases from zero. The crossover from frictionless to 
frictional behavior occurs near /i* « 0.1. We find qual- 
itatively similar behavior for ((f) j) and (z pp ) versus /i e ff 
for the GA model. The GA model gives a 1% larger (cj)j) 
at large /i e ff, which is expected for finite N a . The upper 
right panel of Fig. [2] shows the fraction of floaters Nf /N 
versus \i or /i e ff for the CS and GA models. Both increase 
with /i or /i e ff and then plateau. However, due to slow 
relaxation processes, we detect fewer floaters for the GA 
model, which causes (z pp ) to be 5% larger for the GA 
model at large /i e ff- 

The accumulated mobility distributions (A(() = 

P(x)dx, where ( = |/t|/(/i/ n )) are qualitatively simi- 
lar for the CS and GA models in Fig. 03 At low /i or /i e ff , 
A(Q for both models has a strong peak at £ = 1 [9l.Il3j. 
As ji or /i e ff increases, the peak disappears and the av- 
erage mobility decreases for both models. The quantita- 
tive differences in the mobility distributions for the CS 
and GA models are due to the differences in the tan- 
gential force law shown in the lower left panel of Fig. [TJ 
At fixed amount of overlap, ft/ fn varies linearly with 
u t until the sliding limit at ±u a (and ft/ fn remains 
fixed at ±/i) for the CS model. In contrast, the ratio 

ft/fn = u t /y/(<T?f -rff) 2 -u\ for the GA model is 

periodic in u t with period Au a (as rfj 1 — >• <r^ a ), and thus 
ft/fn decreases when \u t \ exceeds \u a \. Overall, the dif- 
ferences in the structural properties of the MS packings 
obtained from the two models are small. 

However, we identify several distinctive features of the 
structural and mechanical properties of MS packings of 
frictional disks when we consider asperity contacts in the 
GA model. In the lower panel of Fig. [4j we count the 
number of asperity contacts (single, double, and triple) 
for each interparticle contact. We find that MS packings 
from the GA model are isostatic with N£ a = 3N' — 1 
contacts over the entire range of /i e ff- Deviations from 
isostaticity are less than 2% for all N and N a studied. In 
contrast, static packings of frictional particles are hyper- 
static (zpp > 3) when considering interparticle contacts 
for the GA and CS models [4J (c.f lower panel of Fig. [2]). 

By considering asperity contacts, we can also under- 
stand the crossover from frictionless to frictional behav- 
ior in the structural and mechanical properties of static 
packings near /i*. In the top panel of Fig. [4j we plot the 
probability of single and double asperity contacts versus 
/i e ff- We find that single and double asperity contacts 
are roughly equiprobable at low friction, while only dou- 
ble asperity contacts occur at high friction. Thus, to 
maintain isostaticity, at low friction there are typically 
two double and two single asperity contacts per parti- 
cle, while at high friction there are three double asperity 
contacts per particle for a total of approximately six as- 
perity contacts per particle in both cases. We find that 




FIG. 5: The density of vibrational modes D(uj) (in the har- 
monic approximation) for N = 48, N a — 16 and 32, and 
Meff = 1(T 3 , 10~\ 1, and 10 5 for the GA model. The area 
under D(uj) for /x e ff > is the number of nonzero modes 
37V 7 - 2, while for /i eff = it is set to 2N' - 2. The upper-left 
inset shows D(uj) on a log- log scale to highlight low frequen- 
cies. The lower-left inset shows the ratio of the rotational R 
to the translational T content of the modes on a linear fre- 
quency scale. The upper right inset tracks the location cj ma x 
and height i^(cj m ax) of the low- frequency peak in the density 
of vibrational modes versus /i e ff for a system with N a = 32. 
The solid lines have slope —1 and 1. 

the /i e ff where single asperity contacts become less prob- 
able than double asperity contacts (~ 0.1) coincides with 
the characteristic static friction coefficient above which 
the packing fraction, contact number, and mobility dis- 
tributions begin to deviate significantly from frictionless 
behavior. This result has been verified for several N and 
N a . 

The competition between different types of interpar- 
ticle contacts can also be studied in the context of the 
CS model. In the upper panel of Fig. [4j we show the 
probability of low (£ < ( c = 0.5) and high (£ > ( c ) mo- 
bility contacts versus /i. (We find that the results do 
not depend strongly on £ c .) At low friction, most con- 
tacts possess high mobility, while most contacts possess 
low mobility at high friction. In the case of high friction, 
double asperity contacts are analogous to low mobility 
contacts. In the case of low friction, both single and 
double asperity contacts can possess high mobility. As 
in the GA model, the crossover in the probabilities of 
low and high mobility contacts occurs near /i* that sig- 
nals the change from frictionless to frictional behavior in 
((f) j) and (zp P ). 

Another advantage of the GA model is that 
we can directly calculate the density of vibra- 
tional modes D(uo) from the total potential en- 
ergy. The eigenmode with frequency ujj is 

r x,l y,l 9,1 x,N' y,N' 9,N'^ 

rrij {m^ ,ra*' ,m/ , . . . , ra/ ,m|' ,771/ } 

with ^Zxi( m j 1 ) 2 = 1- W e q uan tify the rotational Rj 
and translational Tj content of each mode j, where 



4 



Tj = Ei=i, N > EA=x, y K J ) 2 > and Rj = l- T f , the par- 
ticipation ratio Pj = (Ea^K'T^/^E^K'T) 
for A = x,y and 6 separately, and the optical order 

parameter Q° pt = Y^i,k m^m 9,1 * /(N X^( m j'*) 2 ) tnat 
characterizes whether the rotational content of mode j 
is co- or counter-rotating [lOj. 

The vibrational density of states D(uj) (in the har- 
monic approximation) for frictional MS packings using 
the GA model is shown in Fig. [5j We highlight three 
key features: (i) There is a strong peak at low frequency 
whose height D(uj mdiX ) increases and location co> ma x shifts 
to lower frequency with decreasing /i e ff. We find that 
Wmax ~ Meff and L>(a; max ) ~ as /x eff (cf. upper- 
right inset of Fig. [5]). These modes are mostly rotational 
(R ~ 1), globally incoherent (Q opt ~ 0), and increasingly 
localized (P < 0.1) as /i e ff — >• 0. Similar peaks in D(uo) 
that contain low-frequency rotational modes have been 
found in dimer [HI {N a = 2) and ellipse packings [HI 
at low aspect ratio, which shows that there are impor- 
tant common features between packings of anisotropic 
and frictional particles. For small /i e ff? as uj increases, 



D(uj) approaches the one for frictionless MS disk packings 
with predominantly translational and increasingly local- 
ized modes at high frequencies, (ii) The peak in D(uj) 
at low frequency with R <~ 1 disappears for /i e ff > /i*. 
(iii) For /i e ff > /i*, the modes have mixed rotational 
and translational content with R ~ T at all frequen- 
cies. At low frequencies, the modes are gear-like |20l-[22] 
(Qopt ~ —0.5) and collective (P ~ 0.3). At high frequen- 
cies, the modes are increasingly localized with co-rotating 
angular components (Q Q pt ~ 0-5). 
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